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The shape evolution of two-dimensional islands under electromigration-driven periphery difltusion 
is studied by kinetic Monte Carlo (KMC) simulations and continuum theory. The energetics of the 
KMC model is adapted to the Cu(lOO) surface, and the continuum model is matched to the KMC 
model by a suitably parametrized choice of the orientation-dependent step stiffness and step atom 
mobility. At 700 K shape oscillations predicted by continuum theory are quantitatively verified by 
the KMC simulations, while at 500 K qualitative differences between the two modeling approaches 
are found. 
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I. INTRODUCTION 

The nonequilibrium evolution of two-dimensional 
nanoscale structures is a central theme in current surface 
science, which has been the focus of many experimental 
and theoretical studie o^i^'^'^'^ . On the theoretical side, 
a key challenge is to bridge the length and time scales 
between individual atomic processes and the mesoscopic 
step morphologies that they produce'''^. A highly suc- 
cessful approach to this problem treats the crystal steps 
as continuous entities, whose dynamics is governed by 
the interplay between the adatom population on the ter- 
races, and the processes at the steps which define effec- 
tive boundary conditions on the adatom density fieldi^. 
This approach is well established when the steps are 
rough and sufficiently close to equilibrium for thermody- 
namic considerations to apply at least locally. However, 
atomistic details not captured in the continuum theory 
must generally be expected to become important at low 
temperatures and under strong driving2i2ii2iii. 

In this paper we compare atomistic and continu- 
ous modeling strategies for the specific case of is- 
land electromigrationiSiiiii. The present interest in 
electromigration-induced surface morphology evolution 
is motivated both by its importance for the reliability 
of integrated circuits^^, and by its potential for a rich 
variety of nonlinear dynamic phenomena. Specifically, it 
has recently been shown within a continuum approach 
that the motion and shape evolution of electromigration- 
driven islands in the regime dominated by periphery dif- 
fusion (PD) can be surprisingly complex when the crystal 
anisotropy of the mobility of atoms along the island edge 
is taken into accounljiSiii. Depending on the island size 
and the strength of the anisotropy, a range of different 
dynamic phases involving spontaneous symmetry break- 
ing as well as oscillatory and chaotic shape evolution has 
been predicted. 

With regard to a possible experimental observation of 
these phases, it is important to understand to what ex- 



tent the features of the continuum model survive the dis- 
creteness and stochasticity introduced into the dynam- 
ics on the atomic level. We have therefore conducted 
extensive Kinetic Monte Carlo (KMC) simulations of 
island electromigration, using a realistic model of the 
Cu(IOO) surfacai^. At high temperatures (T = 700 K) we 
find good agreement between the continuum and discrete 
models, and report for the first time (to the best of our 
knowledge) a case of oscillatory shape evolution within a 
fully stochastic KMC simulation. At lower temperatures 
the oscillations persist, but additional features not cov- 
ered by the continuum description appear. The discrete 
and continuum models used in this work are described in 
the next Section. Results at high and low temperature 
are presented in Sect lllll and a concluding discussion can 
be found in Sect II VI 



II. MODELS 

A. Kinetic Monte Carlo model 

Our KMC simulations are based on a lattice-gas 
modeli^ with energetics obtained from effective medium 
theorjii^ (EMT). The energy barriers of the model and 
their relative ordering are in good agreement with ex- 
perimental data for the Cu(OOl) surface^. It should be 
noted that, within the EMT description, barriers on the 
Ag(OOl) and Ni(OOl) surfaces are very similar to those of 
Cu(OOl) up to an overall scaling factor—. The simula- 
tions were implemented using the rejection- free Monte 
Carlo algorithm by Bortz et alm^ with a binary tree 
structure^^. This allowed us to reach time scales on the 
order of seconds in physical time and island sizes up to 
10 000 atoms within reasonable computation time. Our 
model is very similar to that used by Mehl et ali^, who 
also considered the Cu(lOO) surface. However, working at 
somewhat lower temperatures and with smaller islands, 
they did not observe the complex phenomena described 
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in the present paper. 

In principle, the KMC model includes all relevant one- 
particle processes and their corresponding energy barriers 
as given in Ref.l9. However, since our goal is to com- 
pare the KMC simulations with the continuum theory, 
we need to exclude some processes from the KMC simu- 
lations. Since in the continuum model the total mass of 
an island is conserved, we explicitly forbid all adatom de- 
tachment events from an island. This does however not 
exclude a jump into a site without any nearest neighbors 
(NN) but with at least one next nearest neighbor atom. 
This local restriction permits both atoms to jump around 
corner sites and excitation of an atom at the boundary 
with three NN's into an edge atom with one NN through 
two jumps. We also exlude vacancy diffusion inside the 
island, as this is not included in the continuum model. 

It must be emphasized that these restrictions do not 
violate the detailed balance condition in our simulations. 
We have checked that the oscillatory shape evolution phe- 
nomena described in this paper persist even when they 
are relaxed. Needless to say, however, the definition of 
the island size becomes problematic for long times when 
detachment of atoms from the island is allowed for. 

The KMC simulations are conducted as follows. Ini- 
tially an island with a given number of atoms S is placed 
on a lattice in a square or rectangular configuration. The 
hopping rate v of an atom to a vacant NN site on the is- 
land edge is given by 

1/ = 1^0 e-np [- {E s - inm{0, Ann) EB)/kBT] , (1) 

where the attempt frequency fo = 3.0 x 10^^ s^^ and the 
barrier for the jump of an atom at a straight edge is Es — 
0.26 eV. The change in bond number —3 < Ann < 3 is 
the difference in the number of NN bonds between the 
initial and final states with bond energy Eb = 0.26 eV. 
Thus, the model gives an additional kink rounding bar- 
rier E'kr — 0.26 eV (compared to the diffusion along the 
close-packed [110] edges) and a total detachment barrier 
Edct = 0.52 eV from a kink site onto a edge site. 

The electromigration force is introduced by decreasing 
(increasing) the barrier for hops along (against) the di- 
rection of the force by an amount -Ebias- In most of the 
simulations described here the force was directed along 
the a;-direction of the lattice. Experimentally realizable 
values of the biasi^i^ are on the order of E'bias ~ 10~^ 
eV for Cu(lOO) at a current density of 10'' A/cm^. The 
values used in the simulations are typically two orders of 
magnitude larger. However, at least in certain regimes it 
is easy to extrapolate to realistic bias strengths by consid- 
ering larger islands (see the discussion below in Sect lllip . 



B. Continuum model 

The continuum description of two-dimensional 
shape evolution by PD is based on the continuity 
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for the normal velocity w„ of the island boundary. Here s 
denotes the arclength, a the distance between NN sites, 
a the step atom mobility, 7 the edge stiffness, k the local 
curvature, and Ft the tangential component of the elec- 
tromigration force acting on an edge atom. We take the 
force to be of constant strength Fq and aligned along the 
X-axis. With 9 denoting the angle between the normal 
of the island edge and the y-axis (counted positive in the 
clockwise direction), this implies Ft = FqCos{6). 

The crystal anisotropy of the surface enters the con- 
tinuum model through the dependence of 7 and a on 
the edge orientation 9. In previous work based on Q 
the stiffness was assumed to be isotropic, and a simple 
model expression with sixfold symmetry was used for the 
mobilitjiiSiil. Here we are more ambitious, and aim to 
adapt the functions 7(6*) and a (9) as closely as possible 
to the thermodynamics and kinetics of the KMC model. 

The orientation dependence of the step stiffness on 
fcc(lOO) surfaces, in particular on Cu(lOO), has been of 
considerable recent interes1i2i24. For our purposes it is 
sufficient to use the simple expression 
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:, m = tan(0), (3) 



where Ik = {a/2)e'^/^^ is the distance between thermally 
excited kinks, and e — Eb/'^ — 0.13 eV is the kink en- 
ergy. Equation ^ interpolates in a natural way between 
the leading order terms obtained in a low temperature 
expansion of the full expression for the two-dimensional 
Ising lattice gas for m > and m = 0, respectively. For 
the close-packed orientation m = it reduces to the well- 
known relationii^ 
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while for m 7^ and T — s- one obtains the purely en- 
tropic low-temperature expressions^ 
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The transition between the two regimes occurs, and the 
singularity of Q is cut off, when the concentration m 
of forced kinks becomes comparable to the concentration 
a/lk of thermal kinks. 

The anisotropy of the mobility is less well understood, 
and experimental data are so far scarce2S*2&. It is known 
that (7 is isotropic in the absence of an additional kink 
rounding barrienifi*^. Here we extend a phenomenolog- 
ical expression previously derived for the close-packed 
step orientation^'^- to general orientations by writing 
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where pkr = e"^'''/'^^^ is the kink rounding probabiHty. 
The effective kink spacing l^g, which replaces the thermal 
kink spacing Ik at nonzero tilts, is defined in terms of the 
stifi'ness (given by (|2J)) through the relation 



Icb{0) = (aVfcBT)7(0). 



(7) 



This is motivated by the general expression for the rough- 
ness of a step configuration Ci^)^ 



((c(x)-c(x')n = M:| 
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which is valid for steps of arbitrary orientationiiSi^. On 
the basis of © we define the effective kink spacing as that 
length scale at which the amplitude of step fluctuations 
becomes equal to the NN lattice spacing. Equating the 
right hand side of (|HJ| to then immediately leads to 

Q. 

The dependence of the mobility ^ on temperature 
and orientation is fully specified by the values of the 
kink energy and the kink rounding barrier, and repro- 
duces correctly the behavior in various known limits. 
For 9 = the low temperature behavior of the mobil- 
ity depends^ on whether the kink rounding barrier is 
weak (in the sense of i?kr < e) or strong (i?kr > e)- In 



the first case a{0) ~ e 



--Bdot/fcBT 



^ , while in the second case 



For kinked steps les/a is of 
g-(£;dot-i-Bkr)/feBT accor- 



cr(0) ~ e~(^''''*+-^''^"')/''^^ 
order unity, and hence a 
dance with Monte Carlo simulationsiS. For the case of 
a strong kink rounding barrier, which is realized in our 
system, the anisotropy of the mobility [the ratio between 
the maximum and minimum values of cr(0)] is thus of the 
order of e'^l^^^ ^ which is comparable to the anisotropy of 
the stiffness. A detailed comparison of ((BJ with numeri- 
cal measurements of the mobility within the KMC model 
will be presented elsewhere. 

The continuum model ^ can be made dimensionless 
by scaling lengths with 



and time with 
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This implies in particular that increasing the electromi- 
gration force is equivalent (up to a rescaling of time) to 
an increase of the island area by the same factor. Gener- 
ally speaking, electromigration effects become important 
when the linear island size is comparable to Ie- When 
presenting the results for the continuum model the island 
size will be parametrized by the dimensionless radius R 
of a circle of the same area, and time is measured in 
units of tE- For comparison with the KMC simulations, 
the dimensionless radius of an islands of S atoms is then 
defined by ^ = \/sJtt/Ie- 
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FIG. 1: Time sequences of island configurations obtained from 
KMC simulations (left and middle columns) and numerical 
solution of the continuum equation (rightmost column) . Note 
that all configurations are displaced in the y-direction linearly 
in time. Parameters are T = 700 K, S" = 1000, Sbias = 
2 X 10"^ eV (left column), T = 500 K, S = 4000, -Bbias = 10"^ 
eV (middle column), T = 700 K, _R = 2.35 (right column). 



III. RESULTS 
A. 700 K 

In Figure n] we show snapshots of the island shape ob- 
tained from the KMC simulations (left panel at T = 700 
K, middle at T = 500 K) in comparison with the contin- 
uum model results at 700 K (right panel). Apart from the 
shape fluctuations associated with the discreteness of the 
KMC model, the shape evolution at 700 K is seen to be 
very similar to that predicted by the continuum model. 
The dimensionless radius corresponding to the KMC sim- 
ulation parameters is i? = 1.56, somewhat smaller than 
the value i? = 2.35 used in the continuum calculations. 

In the KMC simulations an oscillation starts as ran- 
dom protrusions of the front edge of the island bound- 
ary grow large enough to be amplified by the electromi- 
gration force. The front edge effectively splits into two 
parts, and a straight segment perpendicular to the field 
appears on the island boundary, which is left behind as 
the protrusion advances. In the reference frame of the 
island the standing segment moves backwards and even- 
tually reaches the back end of the island. At this point 
the island has elongated to about twice its initial length. 
When the standing segment disappears the island there- 
fore contracts back into a quasi-stationary shape, from 
which the next oscillation is triggered. The whole se- 
quence of events is identical to that observed in the con- 
tinuum model. The formation of a standing boundary 
segment corresponds precisely to the static facets found 
previously for islands of sixfold symmetrjii^. 

It can be seen from Fig^ that the shape oscillations 
break the symmetry with respect to the horizontal (the 
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FIG. 2: Oscillations of the boundary length obtained from 
KMC simulations (upper panel) and the continuum model 
(lower panel) . The conditions correspond to the left and right 
columns of Fig0 
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FIG. 3: Oscillation period (in units of seconds) measured in 
KMC simulations for islands of size S — 500 — 5000 and bias 
strengths Sbias = 4 x 10"^ eV - 4 x 10"* eV, such that the 
dimensionless island radius is R = 1.56 in all cases. 



direction of the electromigration force) , since the initial 
protrusion may appear either in the upper or the lower 
part of the front edge. The continuum model is strictly 
deterministic. Hence the direction of symmetry break- 
ing is fixed by the initial condition, and each oscillation 
event causes a vertical shift of the center of mass of the 
island in the same direction. This implies that the mean 
direction of island motion forms a nonzero angle with the 
direction of the force, and places this mode of migration 
into the oblique oscillatory (00) phase that was previ- 
ously identified in the case of sixfold symmetrjii^. In 
contrast, in the KMC simulations the direction changes 
randomly from one oscillation event to the next. At least 
for the time sequences of limited length that were gen- 
erated in our simulations, no correlations between sub- 
sequent events could be detected, and hence the island 
center of mass performs a random walk in the transverse 
direction. 

Continuum calculations carried out over a range of di- 
mensionless island radii show that the 00 phase is sand- 
wiched between a narrow regime of oblique stationary 
(OS) motion (where islands of stationary shape move at 
a small angle with respect to the direction of the force) 
for 2.16 < R < 2.24, and a regime where islands break 
up for R > 2.4. For R < 2.16 the motion is stationary 
and straight in the direction of the force (the SS phase). 
Apart from a shift in the locations of the phase bound- 
aries, this is in agreement with the KMC simulations, 
where regimes of stationary and oscillatory behavior fol- 
lowed by island breakup were observed with increasing 
island size. The OS regime could not be unambiguously 
identified in the KMC simulations, possibly because the 
small obliqueness in the direction of motion is easily lost 
in the fiuctuations. 

Quantitatively the oscillation of the island shape can 
be monitored by various observables such as the bound- 
ary length of the island, the average island velocity or the 
center of mass coordinates. In Figure |21 we show the evo- 



lution of the boundary length corresponding to the time 
sequences depicted in Figure ^ Despite a large amount 
of irregularity in the timing of the oscillations observed in 
the KMC simulations, a typical period is clearly visible, 
and the overall shape and amplitude of the individual 
events is quite similar to those obtained from the contin- 
uum model. For the parameters of the KMC simulation 
the characteristic time scale is w 1.3 x 10^* s, which is 
a factor 2.6 larger than the scale used to match the two 
time sequences in Fig|21 Since tE ~ 1%, this discrepancy 
is comparable to that found above in the comparison of 
the effective island radius R. In view in particular of 
the uncertainties associated with the phenomenological 
expression (0 for the step atom mobility, the agreement 
between the two models must be considered remarkably 
good. 

Within the continuum model the behavior of the is- 
land depends only on the dimensionless island size R, 
defined relative to the characteristic length Z^; in © . To 
check if this property carries over to the discrete model, 
we performed a series of KMC simulations in which the 
island size and the strength of the bias were varied in 
such a way that the product S'E^bias remained constant, 
and which therefore correspond to a single value of the 
dimensionless island radius R. As predicted by contin- 
uum theory, the same kind of oscillations depicted in the 
first column of FigQ] was observed for island sizes be- 
tween S = 500 and S = 5000. In addition, the oscillation 
period r was found to scale with island size as a power 
law, 

T^S' (11) 

with z w 1.4 (Fig|21). The measured exponent z is some- 
what smaller than the value z = 2 expected from contin- 
uum theory. There is however a substantial uncertainty 
in the estimates of r both for large islands, where only 
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FIG. 4: Stationary island shapes obtained at T = 700 K with 
the electromigration force directed along the lattice diagonal. 
Full blue lines show time-averaged KMC shapes for island 
size S = 2000 and bias -Ebias = 5 x 10"'' eV (left) and -Ebiaa = 
10~^ eV (right). Dotted red lines show stationary continuum 



shapes of dimensionless radius R 
(right). 
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a small number of oscillations could be observed, and 
for small islands, where the oscillations are obscured by 
noise. 

Exploratory simulations in which the direction of the 
electromigration force was chosen to be different from the 
lattice axes indicate that the oscillatory behavior gener- 
ally persists. An exception is the case when the force 
points along the lattice diagonal, where no oscillations 
are observed. Instead, the islands attain a characteris- 
tic shape which becomes increasingly elongated with in- 
creasing island size or bias strength. This behavior is well 
reproduced by the continuum model. In Fig^the time- 
averaged shapes obtained from the KMC simulations are 
compared to stationary shapes generated by the contin- 
uum model. The match between the two sets of shapes 
is seen to be excellent, provided some adjustement of the 
effective islands size is allowed for. The parameters of 
the KMC simulations in Fig^ imply dimensionless radii 
oi R— 1.106 (left) and R = 1.56 (right), which are again 
somewhat smaller than the values used for the continuum 
model. 



B. 500 K 

As shown in the middle column of Fig^ shape oscilla- 
tions of a type similar to those seen at 700 K are observed 
also in KMC simulations at 500 K. Apart from the re- 
duced amplitude of shape fluctuations, which is due to 
the larger island size, the main difference between the 
left and middle columns is a pronounced elongation of 
the island shape in all phases of the oscillation at 500 K. 

This elongation^^ shows up also in the stationary 
shapes which are characteristic of smaller islands and 
weaker bias (Fig[Sll. A quantitative evaluation of the 
behavior of the aspect ratio of stationary islands as a 
function of island size and bias strength is shown in Fig|Sl 
The aspect ratio is seen to increase with increasing island 



FIG. 5: Time-averaged stationary shapes for KMC islands of 
size S = 500, 1000, 2000 and 3000 atoms at T = 500 K and 
EhisLs = 10~^ eV. For this set of parameters oscillatory shape 
evolution sets in at S = 4000, see the middle column of FiglTI 
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FIG. 6: Island aspect ratio, defined as the ratio of the max- 
imal extensions along and perpendicular to the field, for 
stationary islands obtained from KMC simulations at 500 
K. Left: the aspect ratio as a function of island size with 
Bbias = 7.5 X 10"'' eV, 1.0 x 10"^ eV, and 1.5 x 10"^ eV 
(from bottom to top). The dashed lines are guides to the eye. 
Right: the logarithm of the aspect ratio as a function of bias 
for S = 1000 (lower line) and S = 4000 (upper line). The 
dashed lines are linear fits to the data with the slopes 37 and 
50 for S = 1000 and S = 4000, respectively. 



size up to some upper limit where it saturates. Increas- 
ing the island size further eventually leads to the onset of 
oscillations, as shown in the middle column of FigQ] at 
other parameter values island breakup may also result. 
As was described above, during an oscillation the front 
edge of the island splits up into two parts, one of which 
advances while the other one stays behind. As a conse- 
quence, the island width is halved and the aspect ratio 
doubled at the point of maximal elongation. Between 
oscillations the aspect ratio relaxes back to a value that 
coincides with the saturation value for large stationary 
islands shown in the left panel of FigEI In this sense 
the onset of oscillations can be seen as an attempt of the 
island to increase its aspect ratio beyond the maximal 
value that can be realized by stationary shapes. 

The right panel of FigEl shows that the dependence of 
the island shape on the strength of the electromigration 
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bias is fundamentally different from the dependence on 
island size: Rather than saturating, the aspect ratio of 
stationary islands continues to grow with increasing bias, 
and no onset of oscillations is observed at least for islands 
sizes S < 4000. The inequivalence between size and bias 
dependence implies the breakdown of the continuum the- 
ory, which predicts that the islands shape should depend 
only on the dimensionless parameter S/l"^ ~ S'E'bias- To 
give a quantitative example, the data in Fig|Hlshow that 
an island of size S = 2000 driven by a bias force of 
strength -Ebias = 10^"^ eV attains a stationary aspect ra- 
tio A « 2.2, whereas for S = 1000 and E'bias = 2 x 10"^ 
eV the aspect ratio is twice as large. Within the contin- 
uum theory, the two parameter combinations should lead 
to identical shapes. Indeed, continuum calculations car- 
ried out at 500 K fail to reproduce the elongated shapes 
seen in the KMC simulations. Instead, at the dimension- 
less radius « 1.2 corresponding to the cases described 
above, one finds stationary shapes with an aspect ra- 
tio near unity. Larger islands undergo OS motion for 
2.1 < ^ < 2.2 and break up for R > 2.2. 

Quantitatively the data for the aspect ratio A in the 
right panel of FigElroughly follow an exponential relation 

A - e"-^'"-/'^^'^ (12) 

where a = 37 and 50 for S = 1000 and 4000, respectively. 
While the origin of this behavior so far eludes us, it is 
interesting to note that, in order of magnitude, a ^ \/S. 
This indicates that the relevant energy scale governing 
the behavior of the aspect ratio may be related to the 
energy gain associated with moving a row of atoms (of 
length ~ V^) by one lattice spacing in the direction of 
the electromigration force. 

IV. CONCLUSIONS 

The most important outcome of our study is the ob- 
servation of oscillatory island shape evolution, a phe- 
nomenon first predicted on the basis of continuum 
theorjiiSiii, within a realistic KMC model of the Cu(lOO) 
surface. Since the KMC model has been previously 
validated against experimental observations^, we can 
now confidently pin down the conditions under which 
electromigration-driven shape oscillations should be re- 
alizable in the laboratory. Assuming that a bias strength 
of the order of i?bias = 10~^ eV can be achieved 
experimentalljii^, extrapolation of the KMC results 
shown in Fig|2| indicate that oscillatory behavior of the 
kind shown in Fig^ should emerge at 700 K for islands 
containing about 2 x 10^ atoms, corresponding to an ef- 
fective radius of approximately 90 nm. The period of 
oscillations is expected to lie between 2 s and 40 s, de- 
pending on which value of z is used in the relation Hll|l . 



It is instructive to compare the modes of island mi- 
gration observed in the present study of a surface with 
fourfold crystalline anisotropy to previous work on the 
continuum model, where a sixfold anisotropy of the step 
atom mobility was assumectife. Generally speaking, we 
find that the parameter regime in which complex mo- 
tion occurs (the 00 phase discussed above in Sect lllip 
is less extended, and islands are more prone to breakup 
than in the sixfold symmetric case. Also, the sponta- 
neous symmetry breaking associated with the OS phase 
(which could be detected here only in the deterministic 
continuum dynamics) is more pronounced in the model 
with sixfold anisotropy, in the sense that the angle be- 
tween the direction of motion and the direction of the 
force is much larger—. But it is reassuring that the pre- 
cise form of the anisotropics of stiffness and edge atom 
mobility seems to be inessential for the basic phenomenon 
of oscillatory shape evolution to emerge. 

A detailed explanation of the breakdown of continuum 
scaling and the associated elongation of island shapes ob- 
served in the simulations at 500 K is lacking at present, 
and we can offer only some rather general remarks. In 
previous studies of island shape evolution deviations from 
the continuum approach were typically found when the 
characteristic size of the islands becomes comparable 
to the spacing Ik between thermal kinkgSiiSiii. Since 
electromigration-dominated dynamics sets in at island 
sizes comparable to Ie, it is instructive to consider the 
ratio 

which is seen to increase with decreasing temperature 
and increasing bias strength. For the standard value 
^^bias = 10~^ eV used in our work, the ratio ((T^ in- 
creases from 0.27 at 700 K to 0.49 at 500 K, which makes 
it plausible that the continuum description may become 
questionable in the latter case. According to this cri- 
terion, for an experimentally realizable bias strength of 
-Ebias = 10^^ eV the breakdown would be shifted to 
much lower temperatures, and the continuum descrip- 
tion should remain valid down to room temperature and 
below. 
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